home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / r_correlate.pro < prev    next >
Text File  |  1997-07-08  |  8KB  |  260 lines
  1. ;$Id: r_correlate.pro,v 1.5 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       R_CORRELATE
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes Spearman's (rho) or Kendalls's (tau) rank 
  11. ;       correlation of two n-element vectors. The result is a two-element 
  12. ;       vector containing the rank correlation coefficient and the two-sided 
  13. ;       significance level of its deviation from zero.
  14. ;
  15. ; CATEGORY:
  16. ;       Statistics.
  17. ;
  18. ; CALLING SEQUENCE:
  19. ;       Result = R_correlate(X, Y)
  20. ;
  21. ; INPUTS:
  22. ;       X:    An n-element vector of type integer, float or double.
  23. ;
  24. ;       Y:    An n-element vector of type integer, float or double.
  25. ;
  26. ; KEYWORD PARAMETERS:
  27. ; KENDALL:    If set to a nonzero value, Kendalls's (tau) rank correlation
  28. ;             is computed. By default, Spearman's (rho) rank correlation is
  29. ;             computed.
  30. ;
  31. ;       D:    Use this keyword to specify a named variable which returns the 
  32. ;             sum-squared difference of ranks. If the KENDALL keyword is set
  33. ;             to a nonzero value, this parameter is returned as zero.
  34. ;
  35. ;      ZD:    Use this keyword to specify a named variable which returns the 
  36. ;             number of standard deviations by which D deviates from its null-
  37. ;             hypothesis expected value. If the KENDALL keyword is set to a
  38. ;             nonzero value, this parameter is returned as zero.
  39. ;
  40. ;   PROBD:    Use this keyword to specify a named variable which returns the 
  41. ;             two-sided significance level of ZD. If the KENDALL keyword is 
  42. ;             set to a nonzero value, this parameter is returned as zero.
  43. ;
  44. ; EXAMPLE
  45. ;       Define two n-element vectors of tabulated data.
  46. ;         x = [257, 208, 296, 324, 240, 246, 267, 311, 324, 323, 263, 305, $
  47. ;              270, 260, 251, 275, 288, 242, 304, 267]
  48. ;         y = [201, 56, 185, 221, 165, 161, 182, 239, 278, 243, 197, 271, $
  49. ;              214, 216, 175, 192, 208, 150, 281, 196]
  50. ;
  51. ;       Compute Spearman's (rho) rank correlation of x and y.
  52. ;         result = r_correlate(x, y, d = d, zd = zd, probd = probd)
  53. ;       The result should be the two-element vector:
  54. ;         [0.835967, 4.42899e-06]
  55. ;       The keyword parameters should be returned as:
  56. ;         d = 218.000, zd = -3.64390, probd = 0.000268542      
  57. ;
  58. ;       Compute Kendalls's (tau) rank correlation of x and y.
  59. ;         result = r_correlate(x, y, /kendall)
  60. ;       The result should be the two-element vector:
  61. ;         [0.624347  0.000118732]
  62. ;
  63. ; REFERENCE:
  64. ;       Numerical Recipes, The Art of Scientific Computing (Second Edition)
  65. ;       Cambridge University Press
  66. ;       ISBN 0-521-43108-5
  67. ;
  68. ; MODIFICATION HISTORY:
  69. ;       Written by:  GGS, RSI, Aug 1994
  70. ;                    R_CORRELATE is based on the routines spear.c and kendl1.c
  71. ;                    described in section 14.6 of Numerical Recipes, The Art 
  72. ;                    of Scientific Computing (Second Edition), and is used by 
  73. ;                    permission.
  74. ;-
  75.  
  76. function betacf, a, b, x
  77.   on_error, 2
  78.   eps   = 3.0e-7
  79.   fpmin = 1.0e-30
  80.   maxit = 100
  81.   qab = a + b
  82.   qap = a + 1.0
  83.   qam = a - 1.0
  84.     c = 1.0
  85.     d = 1.0 - qab * x / qap
  86.   if(abs(d) lt fpmin) then d = fpmin
  87.   d = 1.0 / d
  88.   h = d
  89.   for m = 1, maxit do begin
  90.     m2 = 2.0 * m
  91.     aa = m * (b - m) * x / ((qam + m2) * (a + m2))
  92.      d = 1.0 + aa*d
  93.      if(abs(d) lt fpmin) then d = fpmin
  94.      c = 1.0 + aa / c
  95.      if(abs(c) lt fpmin) then c = fpmin
  96.      d = 1.0 / d
  97.      h = h * d * c
  98.      aa = -(a + m) *(qab + m) * x/((a + m2) * (qap + m2))
  99.      d = 1.0 + aa * d
  100.      if(abs(d) lt fpmin) then d = fpmin
  101.      c = 1.0 + aa / c
  102.      if(abs(c) lt fpmin) then c = fpmin
  103.      d = 1.0 / d
  104.      del = d * c
  105.      h = h * del
  106.      if(abs(del - 1.0) lt eps) then return, h
  107.   endfor
  108.   message, 'Failed to converge within given parameters.'
  109. end
  110.  
  111. function gammln, xx
  112.   coff = [76.18009172947146d0,   -86.50532032941677d0,  $
  113.           24.01409824083091d0,    -1.231739572450155d0, $
  114.            0.1208650973866179d-2, -0.5395239384953d-5]
  115.   stp = 2.5066282746310005d0
  116.   x = xx
  117.   y = x
  118.   tmp = x + 5.5d0
  119.   tmp = (x + 0.5d0) * alog(tmp) - tmp
  120.   ser = 1.000000000190015d0
  121.   for j = 0, n_elements(coff)-1 do begin
  122.     y = y + 1.d0
  123.     ser = ser + coff[j] / y
  124.   endfor
  125.   return, tmp + alog(stp * ser / x)
  126. end
  127.  
  128. function ibeta, a, b, x
  129.   on_error, 2
  130.   if (x lt 0 or x gt 1) then message, $
  131.     'x must be in the interval [0, 1].'
  132.   if (x eq 0  or x eq 1) then bt = 0.0 $
  133.   else $
  134.     bt = exp(gammln(a + b) - gammln(a) - gammln(b) + $
  135.              a * alog(x) + b * alog(1.0 - x))
  136.   if(x lt (a + 1.0)/(a + b + 2.0)) then return, $
  137.     bt * betacf(a, b, x) / a $
  138.   else return, 1.0 - bt * betacf(b, a, 1.0-x) / b
  139. end
  140.  
  141. pro crank, w, s
  142.   n = n_elements(w)
  143.   w = [0.0, w]  ;operate on elements w[1], ... , w[n] of the
  144.                 ;n+1 element float array, w.
  145.   s = 0.0
  146.   j = 1
  147.   while(j lt n) do begin
  148.     if(w[j+1] ne w[j]) then begin
  149.       w[j] = j
  150.       j = j+1
  151.     endif else begin
  152.       for jt = j+1, n do begin
  153.         if (w[jt] ne w[j]) then goto, case2
  154.       endfor
  155.       jt = n + 1
  156.       case2:
  157.       rank = 0.5 * (j + jt - 1)
  158.       for ji = j, jt-1 do $
  159.         w[ji] = rank
  160.       t = jt - j
  161.       s = s + t^3 - t
  162.       j = jt
  163.     endelse
  164.   endwhile
  165.   if(j eq n) then w[n] = n
  166.   w = w[1:*]
  167. end
  168.  
  169. function erfcc, x
  170.   z = abs(x)
  171.   t = 1.0 / (1.0 + 0.5 * z)
  172.   erfc = t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+t* $
  173.          (.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+t* $
  174.          (1.48851587+t*(-.82215223+t*.17087277)))))))))
  175.   if x lt 0 then return, 2.0 - erfc $
  176.   else return, erfc
  177. end
  178.  
  179. function r_correlate, x, y, kendall = kendall, d = d, zd = zd, probd = probd
  180.  
  181.   on_error, 2
  182.   nx = n_elements(x)
  183.   if nx le 1 or n_elements(y) le 1 then $
  184.     message, 'x and y must be n-element vectors.'
  185.   if nx ne n_elements(y) then $
  186.     message, 'x and y must be vectors of equal length.'
  187.  
  188.   if keyword_set(kendall) eq 0 then begin ;Spearman's (rho)
  189.     type = size(x)
  190.     wrkx = x
  191.     wrky = y
  192.     ixy  = sort(wrkx) ;Indexes of "wrkx" in ascending order. 
  193.     wrkx = wrkx[ixy]  ;Rearrangement of "wrkx" according to ixy.
  194.     wrky = wrky[ixy]  ;Rearrangement of "wrky" according to ixy.
  195.     crank, wrkx, sf   ;Replace elements of "wrkx" by their rank.
  196.     ixy  = sort(wrky) ;Indexes of "wrky" in ascending order.
  197.     wrkx = wrkx[ixy]  ;Rearrangement of "wrkx" according to ixy.
  198.     wrky = wrky[ixy]  ;Rearrangement of "wrky" according to ixy.
  199.     crank, wrky, sg   ;Replace elements of "wrky" by their rank.
  200.     d = total((wrkx-wrky)^2)
  201.     ;Free intermediate variables.
  202.     wrkx = 0
  203.     wrky = 0
  204.     ixy = 0
  205.     en = nx + 0.0
  206.     en3n = en^3 - en
  207.     aved = en3n/6.0 - (sf + sg)/12.0
  208.     fac = (1.0 - sf/en3n) * (1.0 - sg/en3n)
  209.     vard = ((en - 1.0) * en^2 * (en + 1.0)^2 / 36.0) * fac
  210.     zd = (d - aved) / sqrt(vard)
  211.     probd = erfcc(abs(zd)/1.4142136)
  212.     rs = (1.0 - (6.0/en3n) * (d+(sf+sg)/12.0))/sqrt(fac)
  213.     fac = (1.0 + rs) * (1.0 - rs)
  214.     if (fac gt 0.0) then begin
  215.       t = rs * sqrt((en - 2.0)/fac)
  216.       df = en - 2
  217.       probrs = ibeta(0.5*df, 0.5, df/(df+t^2))
  218.     endif else probrs = 0.0
  219.     ;Return a vector of rank-correlation parameters.
  220.     if type[2] eq 5 then begin 
  221.       d = d+0d & zd = zd+0d & probd = probd+0d
  222.       return, [rs, probrs]
  223.     endif else begin
  224.       d = float(d) & zd = float(zd) & probd = float(probd)
  225.       return, float([rs, probrs])
  226.     endelse
  227.   endif else begin ;Kendall's (tau)
  228.     nnx = 0.0
  229.     nny = 0.0
  230.     is  = 0.0
  231.     ;There seems to be no efficient method of avoiding this nested 
  232.     ;FOR loop structure. An alternate method is possible, but requires 
  233.     ;about (1/2 * nx^2) storage and one FOR loop.
  234.     for j = 0, nx-2 do begin
  235.       for k = j+1, nx-1 do begin
  236.         dx = x[j] - x[k]
  237.         dy = y[j] - y[k]
  238.         aa = dx * dy
  239.         if aa ne 0 then begin
  240.           nnx = nnx + 1
  241.           nny = nny + 1
  242.           if aa gt 0 then is = is + 1 $
  243.             else is = is - 1
  244.         endif else begin
  245.           if dx ne 0 then nnx = nnx + 1 $
  246.           else if dy ne 0 then nny = nny + 1
  247.         endelse
  248.       endfor
  249.     endfor
  250.     d = 0 & zd = 0 & probd = 0 ;Keyword parameters of Spearman's (rho).
  251.     tau = is / sqrt(nnx * nny)
  252.     var = (4.0 * nx + 10.0) / (9.0 * nx * (nx-1.0))
  253.       z = tau / sqrt(var)
  254.     prob = erfcc(abs(z) / 1.4142136)
  255.     ;prob = 1 - errorf(abs(z) / 1.4142136) ;IDL's error function
  256.     return, [tau, prob]
  257.   endelse
  258.  
  259. end
  260.